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We use trapped atomic ions forming a hybrid Coulomb crystal, and exploit its phonons to study an 
isolated quantum system composed of a single spin coupled to an engineered bosonic environment. 
We increase the complexity of the system by adding ions and controlling coherent couplings and, 
thereby, we observe the emergence of thermalization: Time averages of spin observables approach 
microcanonical averages while related fluctuations decay. Our platform features precise control of 
system size, coupling strength, and isolation from the external world to explore the dynamics of 
equilibration and thermalization. 
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How does statistical mechanics emerge from the mi¬ 
croscopic laws of nature? Consider, for example, a finite, 
isolated quantum system: It features a discrete spectrum 
and a quantized phase space, its dynamics are governed 
by the linear Schrddinger equation and, thus, remain re¬ 
versible at all times. Can such a system equilibrate or 
even thermalize? Progress in the theory of nonequilib¬ 
rium dynamics and statistical mechanics sheds light on 
these fundamental questions. It has been shown that 
individual quantum states can exhibit properties of ther¬ 
modynamics depending on entanglement within the sys¬ 
tem [1-7]. While the entire system may very well be 
described by a pure state, any small subsystem and re¬ 
lated local observables may be found in a mixed state 
due to disregarded entanglement with the rest of the iso¬ 
lated system, i.e., the large environment. Further, it is 
predicted that even any individual many-body eigenstate 
of a nonintegrable Hamiltonian yields expectation values 
for few-body observables that are indistinguishable from 
microcanonical averages [8-13]. This conjecture has been 
extensively studied by numerical simulations of specihc 
quantum many-body systems of moderate size, exploiting 
available computational power [14-17]. Recently, there 
have been hrst experiments in the context of thermal¬ 
ization in closed quantum systems with ultracold atoms 
[18-20]. However, fundamental questions on the underly¬ 
ing microscopic dynamics of thermalization and its time 
scales remain unsettled [12, 21, 22]. 

Trapped-ion systems are well suited to study quantum 
dynamics at a fundamental level, featuring unique control 
in preparation, manipulation, and detection of electronic 
and motional degrees of freedom [23-29] . Their Coulomb 
interaction of long range permits tuning from weak to 
strong coupling [30]. Additionally, systems can be scaled 
bottom up to the mesoscopic size of interest to investigate 
many-body physics [31-34]. 

In this Letter, we study linear chains of up to hve 
trapped ions using two different isotopes of magnesium to 
realize a single spin with tunable coupling to a resizable 


bosonic environment. Time averages of spin observables 
become indistinguishable from microcanonical ensemble 
averages and amplitudes of time huctuations decay as the 
effective system size is increased. We observe the emer¬ 
gence of statistical mechanics in a near-perfectly-isolated 
quantum system, despite its seemingly small size. 

The dynamics of our system are governed by the 
Hamiltonian [35, 36] 


Hujz ft-UL 


fkjjjOjaj -\- {a^C -\- a C^) . 

1=1 


( 1 ) 

The spin is described by Pauli operators <Ji {I = x, y, z) 
and h denotes the reduced Planck constant. The hrst 
term can be interpreted as interaction of the spin with 
an effective magnetic held lifting degeneracy of the 
eigenstates of cr^, labeled ||) and |t), while the second 
drives oscillations between these states with spin coupling 
rate Q. The sum represents the environment composed of 
N harmonic oscillators with incommensurate frequencies 
ujj, and the operators aj {oj) annihilate (create) excita¬ 
tions, i.e., phonons, of mode j. The last term describes 
spin-phonon coupling via spin flips, cr^ = {ax ± iay)j2, 
accompanied by motional (de)excitation which is incor¬ 
porated in 


C = exp 


i 


N 


1=1 
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at a strength tunable by Ul and the spin-phonon coupling 
parameters rjj oc IjyluJj. Expanding C in series per¬ 
mits restricting to linear terms for values rjj <C 1 (weak 
coupling). For rjj fa 1, as in our experiment, higher or¬ 
der terms become signihcant (strong coupling), allowing 
the system to explore the full many-body set of highly 
entangled spin-phonon states. This regime is well de¬ 
scribed by full exact diagonalization (ED) only, since 
the discrete nature of the bosonic environment of hnite 
size hinders standard approximations applicable to the 
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spin-boson model considering a continuous spectral den¬ 
sity [36, 47]. 

To study nonequilibrium dynamics of expectation val¬ 
ues (o'i(t)) {I = x,y,z), consider an initial product state 
p{t=0) = p{0) = ps(0) 0 p£;(0), where the spin is in 
a pure excited state, and the bosonic modes are cooled 
close to their motional ground states (average occupa¬ 
tion hj=i,,,Ar < 1). With this, we ensure that energies 
of spin and phonons remain comparable to enable the 
observation of the coherent quantum nature of the dy¬ 
namics which creates entanglement of spin and phonon 
degrees of freedom. Because of the coupling, the spin 
subsystem is in a mixed state for t > 0, even though the 
entire system is evolving unitarily. When thermalization 
occurs, any small subsystem of a large isolated system 
equilibrates towards a thermal state and remains close 
to it for most times [4, 10]. 

The so-called eigenstate thermalization hypothesis 
provides a potential explanation for the emergence of 
thermalization in an isolated quantum system. It can 
be phrased as a statement about matrix elements of few- 
body observables in the eigenstate basis of a many-body 
Hamiltonian [8-13]. Within this conjecture, infinite-time 
averages of expectation values of these observables agree 
with microcanonical averages. A mathematical defini¬ 
tion of this hypothesis and further information are given 
in the Supplemental Material. Based on Refs. [8, 9], to 
interpret experimental results, we assume that a cou¬ 
pling distributes any of the energy eigenstates of an 
uncoupled system {\(j)a)} over a large subset of the 
energy eigenstates of the coupled system {[V'/s)}: i-e., 
\(f)a) = S /3 c,s(o)lV'/ 3 ) [36]. Further, we consider that 
these participating states lie within a narrow energy shell 
around the energy of \4>a) [11, 13]. As introduced in 
Refs. [1, 4, 6, 10], an effective dimension of the subset, 
doff = l/tr(p^), provides an estimation for the ergodicity 
of a system. It has been shown that mean amplitudes of 
time fluctuations of expectation values are bounded by 

1/V^off [4, 6]. 

Correspondingly, for our system, we exploit these pre¬ 
dictions for infinite-time averages, both of spin expecta¬ 
tion values, 

doo((CT/)) = lim - / dt{ai{t)) (3) 

T^OO T Jq 

and of their time fluctuations, 

Soc{{cri)) = \/Poo{{mY) - p,oo{{cri)y- (4) 

To this end, we need to quantify the complexity of the 
dynamics induced by the coupling. Hence, we extend ex¬ 
isting dehnitions of deff to a weighted effective dimension 
[36] 

DeS=^Wa(^\cp{a)\A , (5) 

a ^ P ^ 
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FIG. 1. Complexity of the Hamiltonian studied numerically. 
Parameters are a;i/(27r) = 0.7 MHz, rij=i.,.jv = 1. (a) Di¬ 
mension of truncated Hilbert space dim(77trunc(A^)) for corre¬ 
sponding fractions of initial-state population trptrunc(O) lying 
within "Htrunc (solid lines). For N = S, for example, 85% lie 
within dim('Htrunc) ~ (circle). We derive Deff(A, H, Wj) 
up to dim('Htrunc) = 2^® (dashed line), (b) Choosing Q, and 
varying lJz (dashed line), we can tune the spin-phonon cou¬ 
pling into resonance with different modes (sketched at the 
bottom) and boost the system size. Note, that the actual 
number of participating states is much larger than the normal¬ 
ized quantity Defr; see Eq. (5). (c) For fixed H(A) [cf. dashed 
line in (b)], Defi(a;z) increases significantly with N. This en¬ 
ables the systematic investigation of equilibration and ther¬ 
malization depending on the system size. Error bars show 
systematic numerical uncertainties [36]. 


for p(0) = J2a'^oi\4'a){4‘a\- Here, in contrast to deff, the 
statistical average over Wa is performed after calculating 
the inverse participation ratio for each pure state in the 
mixture [36] . Thereby, Deg also incorporates the number 
of participating states, but is normalized to DeS = 1 for 
the uncoupled system. 

Throughout our Letter, we estimate Dgg numerically. 
DeS depends on N, H, uJz, rji, and p(0). We approx¬ 
imate the latter by truncating the Hilbert space TL to 
^trunc, choosing a phonon number cutoff Uc, such that 
dim('Htrunc) = 2(nc + l)'^ ^ 2^® [36]. For a given compu¬ 
tational power and increasing N, the description of the 
initial-state population by trptrunc(O) becomes less repre¬ 
sentative, leading to increasing systematic uncertainties, 
illustrated in Fig. 1(a). Here, the exponentially growing 
complexity becomes evident: dim('Htrunc) ~ 2^^ is re¬ 
quired to achieve trptrunc(O) = 0.99 for N = 5. Figure 
1(b) highlights the experimental controllability of DeS- 
At {0,012} ~ (2, 1 } X oil, the strong coupling to numer¬ 
ous modes leads to a maximum in DgS- For large o;^, 
the spin can get close to resonance with few modes only, 
the latter composing a comparatively small environment. 
Further, the range of accessible values of DeS grows with 
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FIG. 2. Measured unitary time evolution {a^{t)). Experimental results (black dots, error bars: 1 s.d.) for = 1... 5 compared 
to full ED (solid lines). We exclude numerical results for N = 5 due to their large systematic uncertainties. Oscillations 
(time fluctuations) of high amplitude during the transient duration t/rs < 1, are driven by the evolution of p(0) towards the 
ground state of H. (a) For cuz = 0 and increasing N, excitation is coherently exchanged with a growing number of modes 
resulting in revivals at r^ev (shaded areas), (b) For uj^ ~ Q,{N), expectation values fluctuate around a negative offset. Revivals 
and this nontrivial bias emphasize the coherence of the dynamics, (c) Histograms of experimental measurements sample the 
probability distribution which underlies {o^it)). Here, we show these for t £ [ts, 13rs], Uz ~ 0.{N), and = 1, 5 to exemplify 
the quantities ^exp((crz)) and Se^p{{az)). 


increasing N; see Fig. 1(c). 

We experimentally implement the single spin by two 
electronic hyperfine ground states of ^®Mg+ and add up 
to four ^®Mg+ to engineer the size of the bosonic en¬ 
vironment spanned by N (number of ions) longitudinal 
(axial) motional modes. For details on our experimen¬ 
tal setup, see Refs. [48, 49]. First, we prepare the spin 
state, ps(0) = |i)(4'|, by optical pumping and initial¬ 
ize the phonon state, pe{Q), by resolved sideband cool¬ 
ing [24] close to the ground state. In calibration measure¬ 
ments we determine that the modes are in thermal states 
with ^ 1, which effectively enhances Pj=i...N- 

Next, we apply the spin-phonon interaction by continu¬ 
ously driving Raman transitions with spin coupling rate 
n for variable duration t, where is the controllable 
detuning from resonance [35]. Finally, we detect the spin 
by state-dependent fluorescence. We choose to record 
{az(t))^ while we numerically check that {(Jx,y{t)) feature 
similar behavior. To study dynamics for a large range of 
Heff, we choose 95 parameter settings: We set uji/{2'k) rs 
0.71 MHz which corresponds to an effective spin-phonon 
coupling parameter = rjW^.ni + I « 0.94 for hi = 
1. For each N = 1...5, we use a fixed Q,{N)/(2tt) = 
{0.73(1), 0.95(3), 1.28(3), 1.37(3), 1.58(5)} MHz and vary 
LOz from 0 up to dwi [36]. 

In Fig. 2, we present two sets of {(Jz{t)) ioi N = 1... 5. 
Each data point is obtained by averaging over r = 500 
repetitions yielding an expectation value with statistical 


uncertainty oc If^. We compare {(Tz{t)) with numer¬ 
ical full ED of Eq. (1) with dim(Htrunc) < 2^^. As N 
increases, the accuracy of numerical results decreases sig¬ 
nificantly. For A = 4, we have trptrunc(O) < 0.72. For 
A = 5, since trptrunc(O) < 0.5, we exclude numerical 
results in Figs. 2 and 3; here, even state-of-the-art full 
ED methods [16] could consider trptrunc(O) < 0.75 only 
[36]. For ujz = 0 and A = 1, we confirm oscillations 
of high and persisting amplitude due to the coupling to 
the only mode at wi. For increasing A, the spin cou¬ 
ples to A modes including higher order processes, such 
that spin excitation gets distributed (entangled) into the 
growing bosonic environment. Hence, coherent oscilla¬ 
tions at incommensurate frequencies lose their common 
contrast and appear damped. After the transient du¬ 
ration t/rs ~ 1, with Ts = 27r/H, the spin observable 
remains close to its time average. Still, the conservation 
of coherence of the evolution is evident in our measure¬ 
ments: Revivals of spin excitation due to the finite size of 
the system appear at Trev "^1/ SE, where SE is the mean 
energy difference between modes. And, for ujz ~ nega¬ 
tive time averages of (cTzit)) indicate equilibration of the 
system to the ground state of H, biased by uJz- Both ob¬ 
servations present strong independent evidence that our 
total spin-phonon system is near-perfectly isolated from 
external baths. Independent measurements yield a deco¬ 
herence rate of jdecTs ~ 0.01 [36]. This complements the 
agreement of experimental with numerical results, where 
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Effective magnetic field 


FIG. 3. Time averages and mean amplitudes of time fluctuations of {az{t)). These are calculated from experimental traces (black 
dots, error bars: 1 s.d., derived from the underlying probability distribution of {a^{t)) [36]) for varying and = 1... 5 and 
comparison to full ED (solid lines), (a) Increasing cJz shifts the ground state of H, adjusts spin-mode couplings, and varies Deff. 
Even for small systems, we find agreement of time averages with microcanonical averages, fJ.Bxp{{o'z)) ~ /rmicro((o'z)) (dashed 
lines, shaded areas indicate systematic uncertainties). As D^fi rapidly increases with N, time averages follow microcanonical 
averages for a larger range of lJz. (b) (5exp((o'z)) gradually decreases with N and correlated resonances in /iexp((o'z)) and 
5exp((o'2)) fade away, indicating that we tune our system from microscopic to mesoscopic size. 


we set 7dec = 0. 

We analyze all recorded time evolutions, each contain¬ 
ing S ~ 100 data points in the interval [Ts,13rs], by 
deriving time averages 

Mexp((a4) (6) 

*G[ts,13ts] 


and mean amplitudes of time fluctuations 


<5exp((a’z)) 




tehs.iSTg] 


(7) 

The quantities are illustrated in two examples in 
Fig. 2(c). We plot these in Fig. 3 for fV = 1.. .5 and 
as a function of uJz, together with full ED results for 
N = 1... 4 (solid lines). Tuning ujz across the maximum 
of Deff, cf. Fig. 1(b), and comparing HexpHo'z)) to numer¬ 
ically calculated microcanonical averages Umicm 
(dashed lines) [36] , we find agreement for a larger range of 
LOz when increasing N. This indicates an extended regime 
permitting thermalization. For large lOz, we observe its 
breakdown as the spin couples to an environment of de¬ 
creasing complexity. Finite-size effects are prominent in 
resonances of ^J.expiio'z)) and (5exp((o’z)) for = 1, while 
their amplitudes gradually fade away for higher N. 

For further analysis, we postselect data points well de¬ 
scribed by microcanonical averages, i.e., with a deviation 
of less than 0.1 [36]. For those, we show the dependence 
of SexpUcTz)) on N in Fig. 4(a). Although N sets the size 
of the environment, the complexity of the spin-phonon 
coupling is tuned by D, ojz, rji, p(0), and N, cf. Figs. 1(b) 
and 1(c). Consequently, we study the correlation between 
i5exp((o’z)) and Deff by combining our experimental re¬ 
sults with numerical calculations of DeS in Fig. 4(b). In 


general, mean amplitudes of time fluctuations are pre¬ 
dicted to be upper bounded by 1/V des- For our sys¬ 
tem, we even find a proportionality, i5oo((o’z)) oc 
We motivate this scaling, illustrated by the solid line in 
Fig. 4(b), by a heuristic derivation considering pure ini¬ 
tial states and infinite times, which relies on the eigen¬ 
state thermalization hypothesis [36]. Our measurements 
feature such a scaling for DeS ^ 25, despite our nonide- 
alized initial states and finite observation duration. We 
observe that, for Des > 25, measured mean amplitudes of 
time fluctuations do not further decrease. We attribute 
this to the fact that a system of increasing complexity 
features decreasing energy differences in its spectrum, 
corresponding to smaller relevant frequencies in the dy¬ 
namics. Explicitly, the system requires longer durations 
to approach theoretically predicted values. Here, theory 
considers averages for infinite time, and does not make 
any prediction about relevant time scales in the dynam¬ 
ics. 

In summary, we scale our trapped-ion system includ¬ 
ing its engineered environment up to relevant Hilbert 
space dimensions challenging state-of-the-art full ED. We 
present time averages and fluctuations of a spin observ¬ 
able and exploit an effective dimension to study their de¬ 
pendence on the size of the system. We observe the emer¬ 
gence of quantum statistical mechanics within our iso¬ 
lated system despite its moderate size. Simultaneously, 
we monitor the coherent dynamics of thermalization, re¬ 
vealing the importance of initial and transient time scales 
by direct observation of the evolution towards thermal 
equilibrium. Thereby, we contribute to open questions 
in the field of thermalization [1, 4, 22]. Our approach ad¬ 
mits generating a multitude of initial conditions, choosing 
different system and environment states, and preparing 
initial correlations [24, 25, 27]. In addition, it allows us 
to measure a variety of observables [24, 27, 50] . Applying 
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FIG. 4. Scaling of mean amplitudes of time fluctuations with 
N and Deg. (a) We plot 5exp((nz)) (error bars: 1 s.d.) as 
a function of N. The spread for N < 2 highlights finite-size 
effects, and we show an average value for each N (large sym¬ 
bols, error bars: 1 s.d.). For N = 1 ... we observe a decay 
that ceases for N = 5. (b) 5exp((o'z)) as a function of calcu¬ 
lated Deff (error bars: systematic uncertainties), which cap¬ 
tures the dependence of the effective system size on all exper¬ 
imental parameters. We compare to a scaling Soo oc 1/V Deff 
(solid line), motivated for our system, and our measurements 
agree for Defr < 25. Further increasing Deff, the system needs 
longer durations to resolve decreasing energy differences in 
the environment, unveiling the importance of time scales. 


those techniques, we can study, e.g., non-Markovianity of 
the dynamics, which is evidenced by revivals in the evo¬ 
lution, in detail [51, 52]. Further, increasing the strength 
of the spin-phonon coupling, we can effectively expand 
the observable time span. Possible extensions include 
incorporating more and larger spins, tuning long-range 
interactions, adding external baths [30, 35, 53], and pro¬ 
pelling experimental quantum simulations. Beyond nu¬ 
merical tractability, our experimental setup can be used 
as a test bed to assess the validity of approximated the¬ 
oretical methods that address strong couplings to vibra¬ 
tional baths in a variety of fields, such as molecular and 
chemical physics. 

Recently, we became aware of related studies with 
trapped ions, superconducting qubits, and ultracold 
atoms [54-56]. 
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SUPPLEMENTAL MATERIAL 

TRAPPED-ION HAMILTONIAN 


Our system is described by the Hamiltonian presented 
in Eq. (1), which we recast [35] from: 


H = 


huj. 


N 

i=i 


hn 

'Y 




+ a e 


—iki^z 


)■ 


( 8 ) 

As any two-level system, two electronic states can be 
interpreted as a pseudospin-1/2, represented by the Pauli 
operators ai (/ = x, y, z). Controlling the effective energy 
difference of the states, huJz, via the first term of the 
Hamiltonian is equivalent to applying an effective mag¬ 
netic field to the spin yielding a Zeeman shift that lifts 
the degeneracy of the system in a controlled way. 

The second term describes the environment composed 
of N harmonic oscillators, i.e., bosonic modes, with 
incommensurate frequencies u)j, and their excitations 
(phonons) constitute the environment of the spin in the 
otherwise closed quantum system. In Table I we list de¬ 
scriptions of all parameters in Eq. (8). 

We couple spin and phonons through the photon recoil 
term Here, z is the position of the ion, carrying 

the spin, and fcn represents the effective wave vector of 
the optical field, driving a two-photon stimulated Raman 
transition with spin coupling rate H, the so-called Rabi 
frequency (see below). The related recoil, hk]^, is required 
for creating (annihilating) phonons within the modes 
while flipping the spin, that is, assures energy and mo¬ 
mentum conservation of the spin-phonon coupling. We 
want to emphasize that, in principle, the optical drive in¬ 
duces coherent transitions only, and, consequently, does 
not represent a channel to an external bath. The last 
term in Eq. (8) incorporates both, the term {Ml/2)ax, 
which does not affect the motional state, and the spin- 
phonon coupling + a~C^)^ where the spin- 

flip operators can be related via Ux = (cr’*' -I- a~)j2. 

The position z is expressed as 


N 

kLZ = Vj (a] + Oj) , 


(9) 


j=i 


with corresponding spin-phonon coupling (Lamb-Dicke) 
parameters 


Vj = Mjki, 



( 10 ) 


For values rjj <C 1, we find linear/weak spin-phonon 
coupling that is commonly explored in trapped-ion ex¬ 
periments and 


N 

C'~z^r?j(a]-baj). (11) 

i=i 

We, however, tune r/j ~ 1 to exploit higher-order phonon 
emission/absorption terms and in our experiments the 
full expression. 


C = exp 


N 


3 = 1 


- 1 , 


( 12 ) 


is relevant. Driving at D ~ LOj has two main conse¬ 
quences: First, the coupling of the spin states and the 
spin-phonon coupling get further increased and second, 
the individual transitions between single and multiple 
modes cannot get resolved anymore. Note, that the cou¬ 
pling strength can be effectively enhanced, as rjj is pro¬ 
portional to Aij and, thus, scales with the Fock state 

K): 

Vj ~ + (13) 

In our work, we explicitly make use of the nonlin¬ 
ear/strong coupling and bring our system into a non- 
integrable regime. The system Hamiltonian is exactly 
solvable in the following cases: (i) = 1 and = 0, 

in which it corresponds to the quantum Rabi model [37] , 
and (ii) D = 0, in which spins and phonons are not cou¬ 
pled at all. In any other parameter regime, it can be 
considered a nonintegrable Hamiltonian, by the applica¬ 
ble definitions of nonintegrability [38] : (i) it does not have 
an exact analytical solution, and (ii) there is not a num¬ 
ber of independently conserved observables that equals 
the number of degrees of freedom in the system. 

The full quantum dynamics of the thermalization pro¬ 
cess may be treated with quasi-exact methods, since rel¬ 
evant aspects of thermalization, such as time scales and 
size of fluctuations, critically depend on the set of fre¬ 
quencies involved in the dynamics. In turn, standard ap¬ 
proaches for the spin-boson model, such as path-integral 
methods or generalized master Eq. [39], rely on the as¬ 
sumption that the phonon bath is characterized by a con¬ 
tinuous spectrum of frequencies, and thus, they cannot 
be applied to our mesoscopic system. We access the dy¬ 
namics by full ED, which provides all eigenstates and 
eigenenergies of the total (truncated) Hamiltonian, and 
has been done in recent theoretical studies that investi¬ 
gate similar physics in other systems [11, 14-17, 21]. 


Here, M denotes the mass of the ion chain and Aij is 
the wavefunction amplitude of mode j. Further, it is 
convenient to use rjj = M.j i/wi/wjTyi, with rji describing 
the longitudinal center-of-mass mode (j = 1). 


INVERSE PARTICIPATION RATIO 

The Inverse Participation Ratio (IPR) is used as an 
ergodicity measure to quantify the ability of a system 








parameter name realization physical interpretation variation in experiment 

N number of axial mo- number of ions number of axial harmonic oscillators spanning controlled by deterministic 

tional modes the bosonic environment loading of ions, N = 1... 5 

u)j axial mode frequency mutual Coulomb repulsion resonance frequencies of the harmonic oscilla- controlled by axial confine- 

within axial trapping con- tors (phonon energy hujj) ment, here fixed, m = 

finement 1... N 

Q spin coupling rate Two-Photon Stimulated resonant coupling of the spin states |4.) and If), controlled by laser intensity, 

(Rabi frequency) Raman transition (TPSR) driving coherent oscillations between them chosen for dedicated N to 

maximize Deff 

LOz effective magnetic TPSR detuning from res- lifts degeneracy of coupled spin states |4.} and controlled by variable detun- 
field onant (i.e., carrier) transi- |t), i-e., introduces a bias and alters the total ing 

tion ground state 

r]j spin-phonon coupling momentum kick by the two momentum conservation when creating and controlled via uji (confine- 

(Lamb-Dicke) param- photons of the TPSR annihilating phonons: 77 ^ = (photon recoil en- ment), here fixed, 771 = 0.54 

eters ergy)/(phonon energy); for 77 <C 1 : weak cou¬ 

pling (only carrier transition, i.e., spin flip) 


TABLE I. Detailed description of all parameters in the Hamiltonian. 


to thermalize [4, 10]. Here, thermalization describes a 
process in a closed quantum system, where values of ob¬ 
servables (of any subsystem) equal microcanonical av¬ 
erages [8, 9], while the total system remains in an en¬ 
tangled state during the unitary evolution. In general, 
it requires for (any) few-body observable: (i) the equiv¬ 
alence between time averages and microcanonical aver¬ 
ages, and (ii) decreasing mean amplitudes of time fluctu¬ 
ations around a mean value as a function of the system 
size. 

The IPR is defined for a set of energy eigenstates 
{\ 1 pi 3 )} a-nd an initial pure state as 

IPR(|(/<„)) = ^ I \ (14) 

E/s |c^(a)r 

Recent numerical work shows that systems with increas¬ 
ing values of IPR fulfill condition (i) more accurately, see 
for example Ref. [16], where it is shown that deviations 
of diagonal ensemble predictions from microcanonical en¬ 
semble averages decrease with IPR. Intuitively, we may 
expect that it also leads to condition (ii): We argue that 
a large IPR involves a large set of eigenfrequencies in the 
dynamics of few-body observables and destructive inter¬ 
ferences can lead to a suppression of the size of fluctu¬ 
ations on average. However, large time fluctuations can 
still appear, but during ever shorter durations. Indeed, 
an upper bound on these fluctuations can be obtained 
rigorously based on Ref. [4]. For our system, one can 
find a heuristic derivation that mean amplitudes of time 
fluctuations monotonically decrease as a function of IPR 
(see below). The derivation relies on the Eigenstate Ther¬ 
malization Hypothesis (ETH), plausible assumptions on 
the many-body eigenstates, and a large number of par¬ 
ticipating states (parameter regimes with large IPR). We 
mathematically formulate the ETH and explain it in de¬ 
tail in the last section of this Supplemental Material. To 
examine the validity of our derivation, we perform nu¬ 
merical calculations and compare them to the derived 


scaling: 


Soo'^iicr^)) oc 


1 

IPR(|<('a))’ 


(15) 


which is valid for initial pure states and infinite times. 
In Fig. 5 we show numerical results for 8oo^{{<7z)) for a 
variety of initial conditions and system parameters as a 
function of IPR, and the estimate based on Eq. (15) is 
shown as a solid line. We find reasonable agreement even 
for small values of IPR. 



1/y^lPR(|0„» 


FIG. 5. Numerical study of the scaling of mean amplitudes 
of time fluctuations. Filled and empty symbols correspond 
to initial states \4>o) with one and two phonons per mode, 
respectively. Time fluctuations averaged over infinite time 
are calculated using Eq. (20). Parameters: A'" = 1, = 20, 

fl/(27r) = 0.7 MHz. N = 2, Uc ^ 10, H/(27r) = 1.0 MHz. 
N = 3, ric = 6 , f2/(27r) = 1.3 MHz. N = 4, ric = 5, H/(27r) = 
1.4 MHz. N = 5, Uc = 4, Q,/{2n) = 1.6 MHz. For all values 
oi N, we set a;i/(27r) = 0.7 MHz and take Uz = H/4,..., H/2, 
in steps of H/20. This numerical study substantiates the va¬ 
lidity of our heuristic derivation of Eq. (15) (solid line), and 
our experimental observation of the scaling with effective di¬ 
mension. 

In our study, the initial state is a product state of a 
pure spin state and a motional state cooled close to its 
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ground state, i.e., a mixed phonon state, yet with low 
average occupation. Thus, we need to find a method 
to average the IPR and get an appropriate measure of 
the effective dimension of an initial mixed state. In 
Ref. [6] it is shown that, considering a mixed initial state 
P(0) = '^aWa\(j)a){4>a\, an Upper bound to mean ampli¬ 
tudes of time fluctuations around time-averaged states 
can be found from the definition IPR“^(p(0)) = 
with pp = J2a'^a\cp{o:)\'^. However, this approach has 
the disadvantage that the IPR takes large values for 
highly mixed initial states, even in the case of uncou¬ 
pled systems, since then c/s{a) = and IPR“^(p(0)) = 
choose to define an effective dimension by 
averaging over the IPRs of each state in the mixture, 
which directly leads to our definition in Eq. (5). This 
definition yields Deff = 1 for uncoupled Hamiltonians. It 
takes large values only if interactions lead to a large num¬ 
ber of eigenstates participating in the dynamics, allowing 
us to identify regions where the system should be able to 
thermalize. 

The effective dimension Heff is a weighted measure of 
the number of coupled basis states that are required to 
express the uncoupled basis states. For example, con¬ 
sider a pure state \4>a) in the uncoupled basis, represented 
in the coupled basis {|' 0 / 3 )}, i.e., \4>a) = '^pCp{a)\'il)p). 
Typically, most coefficients |c/ 3 (a)p <C I, such that a sig¬ 
nificant value of Heff is obtained only after considering 
a large number of states each of them contributing 
with a small fraction to the final effective dimension. In 
addition to this, part of our our initial state is a ther¬ 
mal state that incorporates a large number of uncoupled 
states, l^a)- Our definition is chosen such that we do not 
artificially increase Heff by simply adding these states, 
but we renormalize Ueff according to their contribution 
to the mixed state. 

NUMERICAL DIAGONALIZATIONS AND 
THEIR SYSTEMATIC UNCERTAINTIES 

We perform full ED of Eq. (1) to calculate ((T^(t)) and 
corresponding microcanonical averages z)) up 

to Tic- For the calculations, we use Matlab on a work¬ 
station with a quad-core 3.7 GHz processor and 64 GB 
RAM. Besides limitations in computation time, a fun¬ 
damental issue when calculating the full eigensystem of 
a, d X d matrix lies in the size of the available main 
memory that is required to store and process this ma¬ 
trix. Our computer is unable to fulfill this task (al¬ 
ready) for systems of iV = 5 ions with a cutoff at 
ric = 6 phonons per mode that corresponds to a ma¬ 
trix of dimension d^ m 2^® x 2^^. For our experimen¬ 
tal parameters, this forces us to neglect more than half 
of the initial-state population, trptrunc(O) < 0.46, while 
d ~ 1.5 X 2^’^ would be required for trptrunc(O) = 0.75 
and d « 1.4 x 2^^ for trptrunc(O) = 0.99. As a bench¬ 


mark of our numerical effort, we compare this to recent 
theoretical studies that diagonalize maximum sizes of: 
dRi 6.6x10^ Ri 1.0x216 [40], dRi 5.2x10^ Ri 1 . 6 x 2 i 6 [41], 
and d ~ 1.2 X 10^ Ri 1.8 x 2^^ [16], while the authors 
state that those are the maximum sizes tractable by their 
methods. In Table H we list parameters used for all pre¬ 
sented numerical results. 

For calculations of DeS, we use an approximation, i.e., 
a block-diagonalization procedure that exploits the band 
structure of our Hamiltonian matrix. We express Eq. (1) 
in the basis of eigenstates of the uncoupled system \(j)a) = 
ls)lni)i ... jnAr)jv, with s = and \nj) the Fock 

state with n phonons in the mth motional mode. The 
Hamiltonian takes the band matrix form, 

HoL-i,a2 — ( 0 ai|-^| 0 a 2 ) — ,ck 2 4“ k^i,Q 2 ; 

with file matrix of nondiagonal elements. We 

order the states according to their energy such that 
the diagonal entries Ea'^ grow with index number a = 
1... dim('Htrunc)- The matrix Uai.aj couples only those 
states that fulfill —Ea^l ^ \ Va^^a 2 \- Further, we ex¬ 
press the initial state p{t = 0 ) = calcu¬ 

late Dgff of each of the initial states \(j)j) by diagonalizing 

Ha },02 1 and keep only quantum states ai and 0:2 close to 
7 such that (Jai - 7], la 2 - 7I) < fVstates/ 2 , where Wtates 
denotes the number of states kept to calculate the IPR 
of each of the states Thus, H^a \ i® projected onto 
a subspace of Wtates < dim(Htrunc), that optimally con¬ 
tributes to Deff of Finally, we diagonalize Ha}, 02 , 
calculate the contribution of and increase Nstates by 

steps of 1000 , until the resulting values of Des vary by 
less than 1% (5% for N = 5). For example, for Y = 4, we 
use ric = 12 , corresponding to dim('Htrunc) ~ 2 ^®, which 
converges using Ngtates = 5000. However, in the case of 
N = 5 and ric = 7, we are (again) at the limit of our 
computational power, because we need Nstates = 20000 
for converged values of DeS- 

In the following, we describe our procedure to yield 
final values of DeS including a measure of systematic un¬ 
certainties accounting for population which has been ne¬ 
glected by the truncation, trp(O) < 1. We evaluate Dgs 
as a function of truncation 0 < ritcunc < for all param¬ 
eter settings presented in our manuscript (for a list of ric 
see Tab. H). For each parameter setting, we calculate 
the truncated initial-state population trp( 0 , Utrunc), take 
Deff(ntrunc) as a first bound DeS,i, linearly extrapolate 
from it a second bound DeS ,2 with the slope extracted 
from Deff(ntrunc-l), and plot DeS = iDeS,i+DeS, 2)/2 as 
a function of trp(0, ritrunc)- In Fig- 6 , we show an exam¬ 
ple of this. The error bars represent our attributed sys¬ 
tematic uncertainties, given by ± \Dcf{,i — Deff, 2 | /4. For 
N = 1... 4 and all parameter settings, we find that, for 
trp(O) > 0.5, all extracted Deff agree within their uncer¬ 
tainties with corresponding Z)eff(nc), which, in turn, we 
consider most accurate, because for Y < 4, it is ensured 
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N 

ric 

1 

trptrunc(O) 

ric 

2 

trptrunc (0) 

Tic 

3 

trptrunc (0) 

ric 

4 

trptrunc(O) 

ric 

5 

trptrunc{0) 

Fig. 

1(b) 

DeS 





10 

0.94 





1(c) 

DeS 

20 

1.00 

20 

1.00 

20 

1.00 

12 

0.96 

7 

0.72 

2 

(o-z(t)) 

20 

1.00 

20 

1.00 

11 

0.97 

6 

0.72 

(4 

0.35) 

3(a) 

Mmicro((<^2)) 

20 

1.00 

20 

1.00 

16 

0.99 

8 

0.86 

5 

0.46 

3(a) 

Mexp((o'z)) 

20 

1.00 

20 

1.00 

9 

0.93 

5 

0.62 

(3 

0.23) 

3(b) 

<5exp((<^z)) 

20 

1.00 

20 

1.00 

9 

0.93 

5 

0.62 

(3 

0.23) 

4 


20 

1.00 

20 

1.00 

20 

1.00 

12 

0.96 

7 

0.64 


TABLE II. Parameters used for numerical calculations. Cutoff in Fock space ric and truncated initial-state population trptrunc(O) 
used in the numerical calculation for the different Figures. 



FIG. 6. Numerical study of D^s as a function of truncated 
initial-state population. We plot Ileff(u-trunc) for N = 1... 5 
(details see text). The truncation ntrunc = 1 corresponds 
to lowest trptrunc(O), while ntrunc is increased in steps of 
one up to its maximum value ric (limited by our computer 
memory; see Tab. II). The parameters are the same as in 
Fig. 1(c): tJi/(27r) = 0.71 MHz, = 0.54, nj=i...N = 1, 
Q/{2n) = {0.71,0.97,1.21,1.46,1.68} MHz, =_0.8 cui. For 
N = 1.. .4 and for trptrunc(O) > 0.5, we find that Ileff(ntrunc) 
agree within their attributed systematic uncertainties (error 
bars) with numerically most accurate values I?eff(nc) which 
have trp(0, Uc) > 0.96. This suggests that our estimation of 
Deff yields reasonable results, yet with large systematic un¬ 
certainties, even if about 0.5 of the initial-state population 
are neglected in the numerical calculations. 


distribution of the form: 

Mmicro((^2;)) — ^ ^ (^z)/3,/3; (1^) 

P 


where 


P/3 = exp 


{Ep-Ef - 
(AP/2)2 _ 



{Ep-Ef - 
(AP/2)2 _ 


(19) 

The ETH implies that {(Jz)p,p is a smooth function 
of the energy Ep, and, thus, the microcanonical pre¬ 
diction should not depend on the details of the en¬ 
ergy distribution. Further, we apply a procedure sim¬ 
ilar to the one described above for Peff, but based on 
ED (no block-diagonalization), to provide final values of 
Mmicro((u 2 )) including attributed systematic uncertain¬ 
ties. We include numerical calculations of microcanoni¬ 
cal averages in Fig. 3 and Fig. 8, even for N = 5, where 
trptrunc < 0.46. This estimation is only used to detect re¬ 
gions in parameter space where the measured mean value 
is close to the microcanonical average (cp. Fig. 8), and 
is not used in our further analysis. Additionally, numeri¬ 
cal calculations show a faster convergence of ^micro ((^.)) 
with respect to the cutoff ric, since that average is less 
sensitive to the total number of quantum many-body 
states than IPR and Peff- 


that tr/ci(0,nc) > 0.96. This suggests that our procedure 
yields reasonable values for Dgg, yet with large uncer¬ 
tainties, even when up to 0.5 of the initial-state popula¬ 
tion is neglected in the numerical calculation, as it is the 
case for TV = 5. We use for final values Deff = .Deff(nc), 
throughout the manuscript. 

To calculate the microcanonical average, /imicro 
we calculate mean and standard deviation of the energy 
with respect to the initial state, E = (H), and 

AE = \]{H^)-E^. (17) 

We define a microcanonical ensemble by using a Gaussian 


EXPERIMENTAL SETUP 

We employ a linear radio-frequency (rf) Paul trap with 
a drive frequency Di.f/(27r) k, 56 MHz to trap Mg+ iso¬ 
topes with secular frequencies of a;x,y/(27r) « {4.0,4.6} 
MHz (radial direction) and a;i/(27r) k, 0.7 MHz (ax¬ 
ial direction) [48]. Using a photoionization laser (wave¬ 
length A ~ 285 nm), we isotope-selectively load Ix^®Mg+ 
and (TV — l)x^®Mg+, and prepare identical spatial con¬ 
figurations in all measurements [42]. Two electronic 
ground states of the hyperfine manifold of ^®Mg+ (nu¬ 
clear spin / 25 Mg = 5/2) constitute the pseudospin, ]{,) = 
iSi/2\F=3,mF=S) and Jf) = 3S'i/2|T^=2,mF=2), where 
E and mp denote the total angular momentum quan¬ 
tum numbers of the valence electron. Note, that ^®Mg+ 
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f/T-S f/T’S 

FIG. 7. Test for unitarity of the experimentally measured time evolutions, (left) Measured values with statistical uncertainties 
(black dots) and full ED of the unitary time evolution, i.e., with no decoherence (solid line). The shaded area emphasizes 
the time span that is considered to derive /iexp((crz)) and (5exp((crz)) (see Figs. 2 and 3), 1 < f/rs < 13. For t/rs < 15, 
we find our experimental data in good agreement with the numerical calculations, (right) The remaining discrepancy can 
be captured by implying a weak decoherence indicated by the dashed line (a'^{t)) — (az{t))e ydecTs = 0.01, which is 

in agreement with independent calibration measurements. Parameters for the experiments: (a) = 1, {cui, D, a;z}/(27r) = 

{1.92(1), 0.685(5), 0.00(1)} MHz, hi = 0.05(3) (b) N = = (0.73(1), 0.828(5), 0.00(1)} MHz, m = 0.6(1) 

(c) N = 2, (cui, D, cuz}/(27r) = (0.71(1), 0.940(5), 0.00(1)} MHz, n{i, 2 } = (0.4(1), 0.6(1)}. This comparison further justifies to 
describe our system by means of unitary time evolution, and, consequently, treat it as fully isolated within our experimental 
observation. 


has no nuclear spin (/ 26 Mg = 0) and consequently, in our 
realization, it constitutes (only) to the motional mode 
structure. 

To address transitions from 3S'i/2 to either 3Pi/2 or 
3 T 3/2 in ^^Mg+ at wavelengths around 280 nm, we use 
dedicated all solid state laser systems [43, 49]: A laser 
beam with Abd ~ 279.64 nm (labeled BD) is tuned 
rnat/2 below the 35'i/2 to SP 3/2 transition with natural 
line width rp3/2/(27r) = 41.8(4) MHz [44] and aligned 
with a magnetic quantization field |i?| ~ 0.58 inT. The 
cr+-polarized light with intensity Ibd ~ 0.5/sat (satu¬ 
ration intensity Aat — 2500 W/m^) is used for Doppler 
cooling to ~ 1 mK and optical pumping to ](,). This state 
provides efficient, state sensitive detection via the closed 
cycling transition to 3 P 3/2 14,4) to discriminate between 
the ground state manifolds. When BD is shifted near 
resonance (detuning < 27r x 4 MHz), we record photon- 
count histograms with averages corresponding to count 
rates of 100 ms“^ for all F=3 states and 2 ms“^ for 
the F=2 states. We can detect experiments in which 
the ion-chain order is disturbed by a systematic change 
in histogram distributions, and postselect experiments 


where the desired ion-order is achieved, which is true in 
(97.1(3), 82(2), 81(1), 79(2)} % of all our experiments for 
N = 2... 5. Two additional laser beams with Arp « 
280.35 nm (labeled RPl and RP2) are superimposed to 
BD and used for repumping electronic state population 
from 3 Si/2\F=3,mF<3) (RPl) and 35'i/2|F=2) (RP2) 
to II) via 3 Pi/ 2 (with rpi/2/(27r) = 41.3(3) MHz [44]). 
Further, two perpendicular laser beams at ARaman ~ 
279.61 nm (labeled RR and BR) with a /c-vector differ¬ 
ence, fci along the axial direction of the trap (ion-chain 
alignment) enable ground-state cooling of axial motional 
modes. This is achieved via sideband cooling with two- 
photon stimulated Raman transitions (TPSR) detuned 
by A/(27r) = 4-100(10) GHz from the 3S'i/2 to 3 P 3/2 
transition. In the regime of resolved sidebands, e.g., for 
a;i/(27r) = 1.90(1) MHz, TV < 3, and D/(27r) < 0.5 MHz, 
we routinely achieve occupation numbers nj=i...N < 0.1 
by an iterative pulsed sideband cooling procedure. For 
a;i/(27r) = 0.71(1) MHz, we use D/(27r) « 0.25 MHz 
and achieve nj=i,,,N < 1 for TV < 5 by addressing sev¬ 
eral sideband transitions simultaneously. In Table HI, we 
summarize relevant experimental parameters that are de- 
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termined with dedicated calibration measurements, e.g., 
mode-temperature measurements [23]. 


A 

1 

2 

3 

4 

5 

a;i/(27rMHz) 0.724(2) 0.707(2) 0.707(2) 0.708(2) 0.709(2) 

CUjv/oil 

1.00(1) 

1.73(1) 

2.41(1) 

3.05(1) 

3.67(1) 

a;z/(27rMHz) 

[0,1.8] 

[0,2.2] 

[0,2.4] 

[0,2.5] 

[0,2.8] 

H/(27rMHz) 

0.73(1) 

0.95(3) 

1.28(3) 

1.37(3) 

1.58(5) 

A/(27rGHz) 

130(10) 

110(10) 

120(10) 

110(10) 

90(10) 

fii 

0.8(1) 

0.3(1) 

0.6(1) 

0.7(2) 

0.2(1) 

n2 


1.0(2) 

1.1(2) 

1.0(3) 

1.0(2) 

ns 



0.9(2) 

1.0(3) 

1.1(2) 

714 




1.0(4) 

0.7(2) 

ns 





2.3(4) 


TABLE III. Experimental parameters. Lowest (highest) fre¬ 
quency axial vibrational mode effective magnetic 

field ojz varied within interval, spin coupling rate fl, detun¬ 
ing of the Raman beams from the 3Si/2 to 3 P 3/2 transition, 
measured mode occupation number nj=i.,.jv. The detection 
duration is 40 ps for < 2 and 80 /rs for A" > 3. 


All terms incorporating the spin in the Hamiltonian 
are implemented via TPSR as well. In order to achieve 
strong coupling, i.e., Rabi frequencies up to H/(27r) Ri 
1.6 MHz, we use intensities Irr = 2.5(5) x lO^/gat and 
.^BR < 0.8(1) X lO^/gat. We employ two acousto-optical 
modulators to switch on and off beams (turn coupling 
on/off), to tune the relative frequency difference between 
RR and HR, i.e., the detuning of the TPSR (varying Wz), 
and to attenuate the intensity of HR (fine tune H). 

At intensities /rr and /rr, the contribution of sponta¬ 
neous Raman scattering [45] of ^^Mg+ to residual deco¬ 
herence in our experiments, is calculated to yuec < 7(1) 
kHz, assuming a relevant beam waist (1/e^ radius of in¬ 
tensity) w = 55(10) fim for RR and HR. This corre¬ 
sponds to less than 0.07(1) scattering events in t = 13rs, 
resulting in a residual heating effect of less than 0.02 
quanta (in 13rs) distributed among N motional modes. 
This is in agreement with dedicated calibration measure¬ 
ments of the decoherence rate induced by RR and HR. 
In Fig. 7, we compare typical experimental results for 
{(Jz{t)) to corresponding calculated unitary time evolu¬ 
tions with and without considering a decoherence rate 
7dec- From this we infer that TdecTs 0.01 such that 
decoherence contributes to a relative uncertainty of (in 
most cases much less than) 6% to our derived quantities 
during [rs,13rs]. We neglect this effect throughout our 
manuscript and consider our system completely isolated 
from external baths. 



FIG. 8. Comparison of experimental results to microcanonical 
averages. Deviation from /imicro((crz)) (with statistical uncer¬ 
tainties) as a function of Deff (with systematic uncertainties) 
for A = 1... 5. We observe 1 < Defr < 50, while ranges of ex¬ 
perimentally accessible values of for fixed A overlap, and, 
additionally, grow with increasing A. At fixed A, fiez.p{{o'z)) 
approaches /rinicro((cr 2 )) with increasing Deff and the points 
with a deviation < 0.1 (dashed line) are used for evaluation 
in Fig. 4. 


APPLICABILITY OF MICROCANONICAL 
AVERAGES 

In Figure 3(a), we display time averages ^J-exp{{o'z)) 
as a function of effective magnetic field Wz for N = 
1... 5 and compare them with the microcanonical en¬ 
semble averages We correlate the deviation 

lAtexp((o- 2 )) - /rmicro((o- 2 ))| with the corresponding effec¬ 
tive dimension for each data point, see Fig. 8. For 
this sampling of Wz, we find DeS lying between 1 and 
50, e.g., for A = 2 it can be tuned from 2 to 5, while 
for A = 5 the spread is 2 to 50. This implies that we 
can explicitly compare experiments which have equal Heff 
despite distinct values of A. Further, Fig. 8, and ac¬ 
companying numerical calculations, suggest that we can 
continuously scale D^s by increasing A and, particularly, 
tuning Wz. Our definition of D^g measures the ability of 
our system to thermalize, in the sense that fiexpHo'z)) 
approaches ^J,niicYo{{o'z)) as D^g increases. Based on this, 
we select the data points for which we study the corre¬ 
lations between fluctuations and effective dimension in 
Fig. 4: The dashed line in Fig. 8 indicates a deviation 
of 0.1, all measurements closer to /imicro((o' 2 )) are taken 
into account in Fig. 4. 


STATISTICAL UNCERTAINTIES OF TIME 
AVERAGES AND OF MEAN AMPLITUDES OF 
TIME FLUCTUATIONS 

In Figs. 3 and 4, we show time averages, /Xexp((<7z)), 
and mean amplitudes of time fluctuations, (5exp((o’z)), 
both derived from experimentally measured time evo- 
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lutions, ((Tz(<)). In the following, we discuss how we 
obtain statistical uncertainties for both quantities using 
the example of i5exp((o’z))- The quantity i5exp((o’z)), as 
defined in Eq. (7), represents a measure of the mean am¬ 
plitude of time fluctuations. That is, (5exp((<7z)) itself is 
the standard deviation of the expectation values {<7z(t)) 
from their time averaged mean value fJ,exp{{o'z))- The 
statistical uncertainty of (5exp((<7z)), i.e., its standard de¬ 
viation, can be derived from the particular form of the 
underlying probability distribution. To estimate this, we 
use a bootstrapping method: For each experimentally 
measured time evolution (~ 100 data points), we com¬ 
pile 100,000 data sets using the same number of data 
points by sampling with replacement. From these, we 
derive the mean value of ^exp((o'z)) and its standard de¬ 
viation. The latter is then used as statistical uncertainty 
(error bar) for Sexp{{crz))- 

EIGENSTATE THERMALIZATION 
HYPOTHESIS AND MEAN AMPLITUDES OE 
TIME ELUCTUATIONS 

In the following, we give a heuristic derivation of the 
scaling of mean amplitudes of time fluctuations with the 
effective dimension, see Eq. (15), as observed in our ex¬ 
periment. For this purpose, we mathematically formulate 
the ETH and explain it in detail. 

We consider time fluctuations of an operator O. We 
define mean amplitudes of time fluctuations (5oo((0)) (see 
main text), corresponding to the fluctuations of (O) av¬ 
eraged over infinite time. If the system is initially in a 
pure state \(pa) and if it has nondegenerate energy gaps, 
we get [10] 

Soo\{0))= Y. |cftHn%(a)nO;3.Ap, (20) 

01^02 

which depends on both, the structure of the many-body 
eigenfunctions cp{a) and the matrix elements 

We start by discussing the properties of cp{a). Moti¬ 
vated by our experimental conditions, we consider that 
initial states are product states of the form \<j)a) = 
|s)z|?T-i)i • ■ • |?T-Ar) 7 v, where |s)z is a spin state in the 2 - 
basis, and \nj)j is a Fock state with rij phonons in the 
mth vibrational mode. We can write our Hamiltonian as 
H = Hq -\- Hi, where 

hUz X—^ I- t 

Ho = —o-z + 2^ hWja'-aj, 

3 

+ ( 21 ) 

The interaction term Hi (which includes both the spin 
coupling and the spin-phonon coupling term of the spin 
Hamiltonian) couples any given initial state \ 4’a) to differ¬ 
ent states \4>a'), within an energy shell of certain width 


Wq. Thus, we can expect that the many-body eigen¬ 
states, j'i/’/s), will participate in the noninteracting initial 
state only within a range of energies Ep = ± W^. 

Here, and Wa are the mean energy and standard 
deviation of H with respect to jV'a), respectively. This 
observation motivates the ansatz 

cp{a)=FiEp,a)Cp^^, ( 22 ) 

where E{Ep,a) is a smooth function of Ep of width Wa 
defining the energy shell. The function E{Ep, a) satisfies 
the normalization condition ^ D{Ep)\F{Ep,a)\'^dEp = 
1, where we introduce the density of states, D{E) = 
S{E — Ep). Note, D{E) is not representing the effec¬ 
tive dimension, but we will relate the two in the following, 
see Eq. (32). The particular details of Hi will induce a 
factor that is modeled in our ansatz by random variables, 
Ca,p, with average, |C'a,/ 3 P = 1. This ansatz is further 
motivated by random matrix theory: if the interaction 
Hamiltonian, Hi, is a random matrix perturbation, it 
can be shown that the many-body eigenstates satisfy the 
energy shell condition [46] . We expect a similar behavior 
to occur in many-body eigenstates of nonintegrable sys¬ 
tems. Our numerical calculations qualitatively confirm 
this assumption. The wavefunctions appear localized in 
energy-space, as assumed in the energy-shell model, with 
a superimposed random component. We find numerically 
that the energy width takes the value Wa ~ liO for all 
initial states of the form \(j)a) = Yp C/jC'a)]^’/?)- 

The ETH conjecture explains the mechanism behind 
thermalization of nonintegrable closed quantum systems 
by assuming that thermalization is a property of single 
eigenstates. The meaning of this statement can be ex¬ 
pressed as a condition fulfilled by matrix elements of an 
operator O in the eigenstate basis of the nonintegrable 
Hamiltonian [10], 

Opi,p2 = (V’/SijOlV’/Ss) (23) 

= 0{E) 5p^^p^ + -j==f{E,uj)Rp^^p^. 

Here, E = {Ep^ Ep^)l2 and uj = {Ep^ — Ep^), with 
the eigenenergies Ep of our coupled system, and 0[E) 
is a smooth function of E, which corresponds to the mi- 
crocanonical average at energy E. The second term on 
the r.h.s. of Eq. (24) determines the nondiagonal matrix 
elements of O and f{E,uj) is a smooth real function of 
{E, uj). It is centered around w = 0, and it has a typi¬ 
cal width W{E) as a function of w. The set of complex 
numbers Rp^ ^p^ is described by stochastic variables when 
averages are taken over the set of energy eigenstates 

E ^/3i./3J<5/32,/3', (24) 

where the proportionality factor depends on the operator 

O. 
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The normalization of the function f{E,uj) in Eq. (24) 
requires some additional attention. We consider first the 
following summation over matrix elements: 

X! I ./32 I ^ I o I V’/32 ) (#2 IOI #1) 

02 02 

= (#J02|V'/3,) = 1, (25) 

On the other hand, using Eq. (24), we find 

(26) 

02 

^01,02 + X! JJ^\fi^^^)\‘^^*01,02^01<02- 

02 {^0l) ^ ’ 

According to the ETH, is the squared microcanon- 

ical average of operator O for initial states with mean 
energy E. However, focusing on the case in which O is 
a spin operator, ai {I = x,y,z), we note that, in the 
range of parameters in which we observe thermalization, 
we detect small mean values of cr^, and this contribution 
can be neglected. Otherwise, it would lead to a smooth 
dependence of the normalization condition on the vari¬ 
able E, which could be incorporated into the discussion 
below. Using the statistical properties of we get: 

'^\O0^,02?~ X! ^*01,02^01^0-^ 

02 02(¥=0 i) 

« I dE0,D{E0,)-^\f{E,u:)\^. (27) 

Next, we use the approximation D{Ep^) = D{E—uj/2) « 
D{E). This is valid as long as D{E) is a smooth func¬ 
tion that satisfies D"{E)W{EY/D{E) <C 1, such that we 
can substitute D{Eij^) by its average value D{E). Our 
numerical calculations show that this approximation is 
valid within the range of parameters considered in this 
work. Einally, using this approximation together with 
the normalization condition given by Eq. (25), we find 

J dE0,DiEp,)^^\fiE,cu)\^^ J cLo\fiE,uj)\\ 

(28) 

We can also estimate the width of f{E,ijj) in terms of 
the typical width of the energy shell, Wa- For this, let us 
write explicitly the matrix elements of our reference ob¬ 
servable, (Tz, in terms of the unperturbed basis of states, 

I'(’a) I 

(f^^)/3i./32 = C/3 i(Q:i) |(Tz )c^2 (^ 1)5 (29) 

ai 

where we have used explicitly the fact that Cz is diagonal 
in the {(pa) basis. Since C/ 3 j(ai), 0 / 32 ( 01 ) are functions of 
E/ 3 j, Ep^, of width Wai, the energy width of {az)pi,p 2 
will be determined by the typical energy width, Wa^ of 
eigenstates {paP) in the sum. We can even predict from 
Eq. (29) that W{E) « 2Wa « 2nfl. 


Let us use the estimations above in the calculation of 
mean amplitudes of time fluctuations. The latter can be 
expressed, after using the ETH Eq. (24) and the energy 
shell ansatz (Eq. (22)), like 

(5oo^((o-z)) ~ J dEp^dEp^D{EpjD{Ep^) X 

X \F{Ep,,a)p\FiEp,,a)p x 
X ^\f{E,Ep,-Ep,)\\ (30) 

The key observation to evaluate this integral is that both 
F{Ep,a) and f{E, w), have a similar energy width in Ep 
and Lu, respectively, which corresponds to the typical en¬ 
ergy width of the energy shell, Wa- Using the normal¬ 
ization condition for those functions, we can estimate. 


(5oo^((tTz)) OC 


1 


n(-Ea.)Wa 


(31) 


This means that squared time fluctuations approximately 
decay as the inverse of the number of states within the 
energy shell. Assuming, for example, Gaussian distribu¬ 
tions for F{Ep,a) and |/(E,a;)p of width Wa and 2Wa, 
respectively, yields (5oo((o'z)) ~ 0.48/\/D{Ea)Wa- How¬ 
ever, any similar functional dependence that satisfies the 
normalization and energy-shell condition will yield a sim¬ 
ilar scaling with D{Ea)Wa. 

We can finally relate mean amplitudes of time fluctu¬ 
ations to the IPR by noting that, 


IPR(|</>a)) 


1 


1 

lE,D{Ep)\F{Ep,a)P 
CX D{Ea)Wa. 


(32) 


This relation is, again, a result of the normalization and 
energy shell condition of the function F{Ep, a). It has a 
very physical interpretation, since it relates the IPR to 
the number of states within the energy shell defined by 
Wa- Combining the last two equations, we arrive at 


docP^{{(Tz)) OC 


1 

IPR(l<('a))' 


(33) 


Our numerical calculations confirm the scaling predicted 
by Eq. (33). In Fig. 5 we calculate numerically the mean 
amplitudes of time fluctuations by using Eq. (20) for a 
variety of initial conditions and system parameters. Our 
heuristic derivation assumes that the system is ergodic 
enough, so that a large number of many-body eigenstates 
participate in the initial state. Thus, to check Eq. (33), 
we choose a wide range of parameters where the system 
is observed to thermalize efficiently. Figure 5 shows that 
the estimate given by Eq. (33) works surprisingly well, 
even for small numbers of particles. 












